
;+
;PURPOSE
;	Code that reads the cluster file from SLUG
;SYNTAX
;	read_slug, [modelname, outputname=outputname, /here, /slow, path=path,
;			/skip_cluster,
;			hist_outputname=hist_outputname, int_outputname=int_outputname]
;INPUTS
;	modelname: the name of the model as determined in your slug.param
;		   file when you ran the slug code. This code requires you
;		   to have the $SLUG_DIR environment variable.
;		   note that the default value for slug.param is SLUG_DEF
;	outputname: name of output file [default=modelname+'.fits']
;		this is the name of a fits file that gets written
;		(extension 1) is the array of structures for each cluster
;			mass: mass of cluster in solar masses. this is the
;				sum of the masses of the stars
;			age: the age of the cluster when evaluated
;			time: the time from the start of the simulation
; 	     additional tags: these are the fluxes named by the names of the
;				fluxes used
;	hist_outputname: name of output fits file that parsed histograms 
;				[default='histo_'+outputname]
;	int_outputname: name of output fits file that parsed integral photometry
;				[default='integral_'+outputname]
;	imf_outputname: name of output fits file that parsed imf historam file
;				[default='imf_'+outputname]
;	misc_outputname: name of output fits file that parsed misc. file
;				[default='misc_'+outputname]
;      path:     the location of the model. This has precendence over /here
;KEYWORDS
;	/here: if set then read_slug will read and write the 
;		file from the current folder
;	/slow: set to run slower but with better memory management
;	/skip_cluster:set to skip the parsing of the cluster file
;	/imf_skip: set to skip parsing IMF file
;	/misc_skip: set to skip parsing misc file
;NOTES
;	note that the seed is written to the header
;
;NOTE!!!: the field galaxy bin is stored in a cluster with mass of -1 and
;	has one entry for each timestep
;
;USEFUL NOTE FOR DEALING WITH THIS PROGRAM'S OUTPUT
;	if the output file is large you can read in certain elements with 
;	mrdfits() through use of the rows keyword
;
;Written by R. da Silva, UCSC, 4-21-2010
;Modified '/here' by MF, Feb 2011 
;Added '/slow' keyword, RdS, Feb 2011
;-

pro read_slug, filename, outputname=output_name,int_outputname=int_outputname,$
         misc_outputname=misc_outputname,$
	 hist_outputname=hist_outputname, here=here, path=path,slow=slow,$
	skip_cluster=skip_cluster, imf_outputname=imf_outputname, $
	imf_skip=imf_skip, misc_skip=misc_skip
;-------------------------------------------------
;	Initialize some stuff
;-------------------------------------------------

;check inputs
if NOT keyword_set(filename) then begin
	print, '***Filename is not set***'
	print, 'Using default = SLUG_DEF'
	filename='SLUG_DEF'
endif
;set defaults

if keyword_set(here) and not keyword_set(path) then path='./' 
if not keyword_set(here) and not keyword_set(path) then path='$SLUG_DIR/output/'
if not keyword_set(output_name) then output_name=path+filename+'.fits'

;;;if keyword_set(skip_cluster) then slow=1
;check if 64 or 32 bit machine since it makes a difference how your
;longs are allocated and their datatypes
machine_64bit=is64bit() 
if machine_64bit then longdt=14 else longdt=3
hdr=strarr(5)


;-------------------------------------------------
;	Read in the guide file so I know how to  |
;	parse the file				 |
;-------------------------------------------------

;if keyword_set(here) EQ 0 then openr, lun1,'$SLUG_DIR/output/guide_'+filename+'.dat', /get_lun else $
;	openr, lun1,'./guide_'+filename+'.dat', /get_lun
if x_chkfil(path+'guide_'+filename+'.dat', /silent) then begin
openr, lun1, path+'guide_'+filename+'.dat', /get_lun
checkingit=fstat(lun1)
if checkingit.size GT 0 then begin


;-----------Begin reading in data to put in hdr----------
;read in the which_fluxes bit mask
which_fluxes=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]

;read in the Seed (long)
seed=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
sxaddpar, hdr, 'SEED', seed, $
	'seed used for simulation'
;read in ntime (long)
max_time=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
sxaddpar, hdr, 'MAXTIME', max_time, $
        'the model was run until this time [yr]'
;read in timestep (float)
timestep=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
sxaddpar, hdr, 'TIMESTEP', timestep, $
        'timestep between evaluations [yr]'
;read in galaxy_sfr (float)
sfr=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
sxaddpar, hdr, 'SFR', sfr, $
        'star formation rate of galaxy [Msun/yr]'
;read in imf_type (long)
imf=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
case imf of
 0:imf_name='Kroupa'
 1:imf_name='Salpeter'
 2:imf_name='Chabrier'
 3:imf_name='IGIMF'
 4:imf_name='SLOPE'
endcase
sxaddpar, hdr, 'IMF', imf_name, $
        'stellar initial mass function used'
;read in metallicity (float)
metal=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
sxaddpar, hdr, 'METAL', metal, $
        'Z: note Z_solar=.020'

;;make a master header
hdr_mast=hdr

;-------------end read in data for header-----------

gfstats=fstat(lun1)
;find out size of file and subtract out stuff I've already
;read
gfstatsize=gfstats.size-n_bytes(which_fluxes)-n_bytes(seed)-n_bytes(max_time)-$
	n_bytes(timestep)-n_bytes(sfr)-n_bytes(imf)-n_bytes(metal)

gndoubles=gfstatsize/8;number of doubles remaining in the whole file
if keyword_set(slow) EQ 0 AND gndoubles GT 0 then begin

  arr=read_binary(lun1,DATA_TYPE=5,DATA_DIMS=gndoubles)
  wh=dindgen(gndoubles/3d)*3d
  mass=arr[wh] ;mass of every cluster
  how_many=arr[wh+1] ;how many timesteps it was evaluated at
  first_time=arr[wh+2] ;what timestep was its first
  id=lindgen(n_elements(mass))
  ;delvarx, wh
  undefine,wh
endif
close, lun1
free_lun, lun1
;-------------------------------------------------
;       Read in the output file		 	 |
;-------------------------------------------------

;open file
;if keyword_set(here) EQ 0 then openr, lun, '$SLUG_DIR/output/'+filename+'.dat', /GET_LUN else $
;	openr, lun, './'+filename+'.dat', /GET_LUN 

openr, lun, path+filename+'.dat', /GET_LUN 

;-----------------------------
;read in number of time steps
;-----------------------------
ntime=(read_binary(lun,DATA_TYPE=longdt,DATA_DIMS=1))[0]
fstats=fstat(lun)
fstatsize=fstats.size-n_bytes(ntime)
ndoubles=fstatsize/8;number of doubles remaining the whole file
if keyword_set(slow) EQ 0 then begin
  arr=read_binary(lun,DATA_TYPE=5,DATA_DIMS=ndoubles)
  close, lun
endif

;----------------------------------------
; figure out which tables are written in
;----------------------------------------
all_tables=['NUV', 'FUV', 'BOL', 'QH0',$
'J','H','K','SDSSu'$
,'SDSSg','SDSSr','SDSSi'$
,'SDSSz','johnsonU','johnsonB'$
,'johnsonV','johnsonR'$
,'johnsonI']
ntables=0L
for wh_tab=1, n_elements(all_tables) do begin
  if ((which_fluxes AND 2L^wh_tab)  NE 0) then begin
     if ntables EQ 0 then tab_ind=wh_tab-1 else tab_ind=[tab_ind, wh_tab-1]
     ntables++
  endif
 endfor


;---------------
;now I need to parse with where statements
;---------------
if keyword_set(slow) then begin
	        arr=read_binary(lun,DATA_TYPE=5,DATA_DIMS=ntime) 
   	        time_steps=arr ;the time array 
	endif else begin
		time_steps=arr[0:ntime-1] ;the time array 
		if n_elements(arr) EQ ntime then skip_cluster=1
		if keyword_set(skip_cluster) EQ 0 then arr=arr[ntime:n_elements(arr)-1]
	endelse

if keyword_set(skip_cluster) EQ 0 then begin ;;; I am moving this up
if keyword_set(slow) EQ 0 then ncluster=total(how_many)


tagnames=['id',$
'mass', $
'mass_actual', $
'age_clus', $
'age_gal', $
'nstars', $
'max_stellar_mass']
n_nonflux_tags=n_elements(tagnames)
tagnames=[tagnames, all_tables[tab_ind]]
descrip=['D(1)', replicate('D(1)', 6+ntables)]

rcreate_struct, str, 'SLUG_'+filename, tagnames, descrip

if keyword_set(slow) EQ 0 then str=replicate(str, ncluster)
;below is a complicated series of steps required to get the right values in the right places
if keyword_set(slow) then begin
  openr, lun1, path+'guide_'+filename+'.dat', /get_lun
  dum=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
  dum=(read_binary(lun1,DATA_TYPE=4,DATA_DIMS=1))[0]

  gfstats=fstat(lun1)
  ;find out size of file and subtract out stuff I've already read
  gfstatsize=gfstats.size-n_bytes(which_fluxes)-n_bytes(seed)-n_bytes(max_time)-$
        n_bytes(timestep)-n_bytes(sfr)-n_bytes(imf)-n_bytes(metal)
  gndoubles=gfstatsize/8;number of doubles remaining in the whole file
endif


;;if keyword_set(skip_cluster) EQ 0 then begin ;;; I am moving this up

age_galaxy=0d
ss=systime(3)
ipart=0
writ=0
done=0L
entry=0L
if keyword_set(slow) then niter=gndoubles/3d else niter=n_elements(mass)
for i=0L, niter-1 do begin
  if keyword_set(slow) then begin
     arr=read_binary(lun1,DATA_TYPE=5,DATA_DIMS=3)
     mm=arr[0]
     hm=arr[1]
     ft=arr[2]
     id_=i
     done=0L
     entry=0L
     arr=read_binary(lun,DATA_TYPE=5,DATA_DIMS=long(4+ntables*hm))  
  endif else begin
	mm=mass[i]
	hm=how_many[i]
	ft=first_time[i]
	id_=id[i]
  endelse
  if (i mod 5d5 EQ 0 AND i NE 0) OR i EQ n_elements(mass)-1 then begin
	print, 100*float(i)/niter, $
		' percent completed', memory(/current)/1d6, (systime(3)-ss)/60.
  endif
  m1=arr[entry];read mass
  entry++ 	;increment entry
 nstars=arr[entry]  ;;
 entry++;;
 biggest_star=arr[entry];;
 entry++;;
 age_galaxy=arr[entry]
 entry++
  if hm EQ 0 then continue
  if keyword_set(slow) then  str=replicate(str[0], hm)
  for j=0., hm-1. do begin
	str[done].nstars=nstars;;
	str[done].max_stellar_mass=biggest_star;;
	str[done].id=id_
 	str[done].mass_actual=m1
	str[done].mass=mm
	str[done].age_gal=time_steps[ft+j]
	str[done].age_clus=str[done].age_gal-age_galaxy

	for ntab=0, ntables-1 do begin
	  str[done].(n_nonflux_tags+ntab)=arr[entry]
	  entry++
	endfor
	
  done++
endfor

if keyword_set(slow) then begin
  if writ EQ 0 then begin
	writ=1
	mwrfits, str, output_name, hdr, /create 
	endif else $
	mwrfits_chunks, str, output_name, /append, /silent
endif
endfor
endif
if keyword_set(slow) then begin
close, lun
free_lun, lun
close, lun1
free_lun, lun1
endif else mwrfits, str, output_name, hdr, /create
if keyword_set(arr) then undefine, arr
;-----------------------------------------------------------------------------------
;		PARSE histogram file						   |
;-----------------------------------------------------------------------------------
if not keyword_set(hist_outputname) then $
	;hist_outputname='$SLUG_DIR/output/histo_'+filename+'.fits'
         hist_outputname=path+'histo_'+filename+'.fits'

;if keyword_set(here) EQ 0 then $
;	 openr, lun1,'$SLUG_DIR/output/histo_'+filename+'.dat', /get_lun else $
;        openr, lun1,'./histo_'+filename+'.dat', /get_lun 
openr, lun1, path+'histo_'+filename+'.dat', /get_lun 


nx=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
if nx NE 0 then begin

ny=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
ntime=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
nhist=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]

if neareq(time_steps[0], 0) then tm=-1 else tm=0

histo=dblarr(nx,ny,ntime+tm)
str={histo:histo,$
log_x_axis:dblarr(nx),$
log_y_axis:dblarr(ny),$
;mag_x_axis:dblarr(nx),$
;mag_y_axis:dblarr(ny),$
x_axis_name:'',$
y_axis_name:'',$
time_axis:dblarr(ntime+tm)}
for hwc=0,nhist-1 do begin

hf_ind1=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
hf_ind2=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]

str.y_axis_name=(all_tables[tab_ind])[hf_ind2]
str.x_axis_name=(all_tables[tab_ind])[hf_ind1]+'-'+(all_tables[tab_ind])[hf_ind2]

xmin=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]
xmax=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]
xbin=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]

ymin=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]
ymax=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]
ybin=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=1))[0]

str.log_x_axis=mkarr(xmin, xmax-xbin, xbin)
str.log_y_axis=mkarr(ymin, ymax-ybin, ybin)

;str.mag_x_axis=-2.5*(str.log_x_axis)-48.6
;str.mag_y_axis=-2.5*(str.log_y_axis)-48.6


histo1=(read_binary(lun1,DATA_TYPE=5,DATA_DIMS=nx*ny*ntime))
histo1=reform(histo1, ntime,ny, nx)
histo1=transpose(histo1,[2,1,0])

if neareq(time_steps[0], 0) then begin
	str.time_axis=time_steps[1:*] 
	histo1=histo1[*,*, 1:*] 
	endif else str.time_axis=time_steps ;this else should never be reached for current
						;implementation of the code
str.histo=histo1
if hwc EQ 0 then str1=str else str1=[str1, str]
endfor

hdr=hdr_mast
mwrfits, str1, hist_outputname, hdr, /create

endif
close, lun1
free_lun, lun1
;-----------------------------------------------------------------------------------
;		PARSE integral file						   |
;-----------------------------------------------------------------------------------
if not keyword_set(int_outputname) then $
         int_outputname=path+'integral_'+filename+'.fits'
openr, lun2, path+'integral_'+filename+'.dat', /get_lun
ntime=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]
ntables=(read_binary(lun1,DATA_TYPE=longdt,DATA_DIMS=1))[0]

tagnames2=['age_gal']
n_nonflux_tags2=n_elements(tagnames2)
tagnames2=[tagnames2, all_tables[tab_ind]]
descrip2=[replicate('D', 1+ntables)]

rcreate_struct, str2, 'SLUG_INT_'+filename, tagnames2, descrip2
str2=replicate(str2, ntime)
kk=0
fstats2=fstat(lun2)
fstatsize2=fstats2.size-n_bytes(ntime)*2
ndoubles2=fstatsize2/8;number of doubles remaining the whole file
arr2=read_binary(lun1,DATA_TYPE=5,DATA_DIMS=ndoubles2)
ee=0
for tc=0, ntime-1 do begin
  str2[tc].age_gal= time_steps[tc]
  for fc=1, ntables do begin
    str2[tc].(fc)=arr2[ee]
    ee++
  endfor
endfor
close, lun2
free_lun, lun2
if n_elements(str2) GT 1 and neareq(str2[0].age_gal, 0) then str2=str2[1:*]
hdr=hdr_mast
mwrfits, str2, int_outputname, hdr, /create

;-----------------------------------------------------------------------------------
;                    PARSE IMF file                                                |
;-----------------------------------------------------------------------------------
if not keyword_set(imf_skip) then begin

if not keyword_set(imf_outputname) then $
         imf_outputname=path+'imf_'+filename+'.fits'
openr, lun3, path+'imf_'+filename+'.dat', /get_lun
if EOF(lun3) EQ 0 then begin ;check if file is empty
nbins=(long(read_binary(lun3,DATA_TYPE=5,DATA_DIMS=1)))[0]

bins=(read_binary(lun3,DATA_TYPE=5,DATA_DIMS=nbins))

tagnames3=['id', 'imf_hist']
descrip3=['D', 'D('+rstring(nbins)+')']

rcreate_struct, str3, 'SLUG_IMF_'+filename, tagnames3, descrip3
fstats3=fstat(lun3)
fstatsize3=fstats3.size-n_bytes(bins[0])*(nbins+1)
ndoubles3=fstatsize3/8;number of doubles remaining the whole file
nobj=ndoubles3/nbins
hdr=strarr(nbins+1)
for i=0, nbins-1 do sxaddpar, hdr, 'BIN'+rstring(i), bins[i], $
	'Right side of bin for IMF_HIST in solar masses'

writ=0
if keyword_set(slow) then begin
  for i=0L, nobj-1 do begin
     arr=(read_binary(lun3,DATA_TYPE=longdt,DATA_DIMS=nbins))
     str3.id=i
     str3.imf_hist=arr
     if keyword_set(writ) EQ 0 then begin
	writ=1
	mwrfits, str3, imf_outputname, hdr, /create 
     endif else mwrfits_chunks, str3, imf_outputname, /append
  endfor
endif else begin
  str3=replicate(str3, nobj)
  arr=(read_binary(lun3,DATA_TYPE=5,DATA_DIMS=ndoubles3))
  str3.id=lindgen(nobj)
  str3.imf_hist=reform(arr, nbins, nobj)
  mwrfits, str3, imf_outputname, hdr, /create
endelse

mwrfits, bins, imf_outputname
endif;close check if file is empty
close, lun3
free_lun, lun3

endif

;-----------------------------------------------------------------------------------
;                    PARSE MISC file                                               |
;-----------------------------------------------------------------------------------
;currently done even if guide file is empty
if keyword_set(misc_skip) EQ 0 then begin
if not keyword_set(misc_outputname) then $
         misc_outputname=path+'misc_'+filename+'.fits'
openr, lun4, path+'misc_'+filename+'.dat', /get_lun
arr=(read_binary(lun4,DATA_TYPE=5,DATA_DIMS=24))
sub=lindgen(23)
hist_sub=where(((sub-1) mod 2) EQ 0)
misc_str={total_galaxy_mass:arr[0], cmf_hist_mass:arr[hist_sub],$
	  cmf_hist_num:arr[hist_sub+1], sfh_hist_t:dblarr(arr[23]),$
	sfh_hist_mass:dblarr(arr[23])}
misc_str.sfh_hist_t=1d6*findgen(arr[23])
arr=(read_binary(lun4,DATA_TYPE=5,DATA_DIMS=long(arr[23])))
misc_str.sfh_hist_mass=arr

close, lun4
free_lun, lun4
if keyword_set(hdr_mast) EQ 0 then hdr_mast=strarr(1)

mwrfits, misc_str, misc_outputname, hdr_mast, /create
endif
endif else print, 'Guide file has zero size... no star born...'
endif else print, 'Guide file empty... no stars born...'

;;just in case free
if keyword_set(lun) then free_lun, lun
if keyword_set(lun1) then free_lun, lun1
if keyword_set(lun2) then free_lun, lun2
if keyword_set(lun3) then free_lun, lun3
if keyword_set(lun4) then free_lun, lun4

end
